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ABSTRACT 

Many astrophysical applications involve magnetized turbulent flows with shock waves. Ab initio star for- 
mation simulations require a robust representation of supersonic turbulence in molecular clouds on a wide 
range of scales imposing stringent demands on the quality of numerical algorithms. We employ simulations 
of supersonic super-Alfvenic turbulence decay as a benchmark test problem to assess and compare the per- 
formance of nine popular astrophysical MHD methods actively used to model star formation. The set of nine 
codes includes: ENZO, FLASH, KT-MHD, LL-MHD, PLUTO, PPML, RAMSES, STAGGER, and ZEUS. 
These applications employ a variety of numerical approaches, including both split and unsplit, finite difference 
and finite volume, divergence preserving and divergence cleaning, a variety of Riemann solvers, a range of 
spatial reconstruction and time integration techniques. We present a comprehensive set of statistical measures 
designed to quantify the effects of numerical dissipation in these MHD solvers. We compare power spectra for 
basic fields to determine the effective spectral bandwidth of the methods and rank them based on their relative 
effective Reynolds numbers. We also compare numerical dissipation for solenoidal and dilatational velocity 
components to check for possible impacts of the numerics on small-scale density statistics. Finally, we dis- 
cuss convergence of various characteristics for the turbulence decay test and impacts of various components of 
numerical schemes on the accuracy of solutions. The nine codes gave qualitatively the same results, implying 
that they are all performing reasonably well and are useful for scientific applications. We show that the best 
performing codes employ a consistently high order of accuracy for spatial reconstruction of the evolved fields, 
transverse gradient interpolation, conservation law update step, and Lorentz force computation. The best re- 
sults are achieved with divergence-free evolution of the magnetic field using the constrained transport method, 
and using little to no explicit artificial viscosity. Codes which fall short in one or more of these areas are still 
useful, but they must compensate higher numerical dissipation with higher numerical resolution. This paper 
is the largest, most comprehensive MHD code comparison on an application-like test problem to date. We 
hope this work will help developers improve their numerical algorithms while helping users to make informed 
choices in picking optimal applications for their specific astrophysical problems. 

Subject headings: ISM: structure — Magnetohydrodynamics: MHD — methods: numerical — turbulence 
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L EvlTRODUCTION 

It is well established that the observed supersonic turbu- 
lence plays an important role in the fra gmentation of molec- 
ular clouds leading to star for mation (iMac Low & KlessenI 
I2004t iMcKee & Ostrikej|2007h . As illustrated by numeri- 
cal simulations, random supersonic flows in an isothermal 
gas result in a complex network of shocks creating a fil- 
amentary density structure with a very large dens ity con- 
trast (e.g.. iKritsuk et al.ll2007HFederradi et aDl2008l see also 
references in iKlessen et all (l2009h : iPudritzl d201lb ). Be- 
cause it can naturally generate density enhancements of suf- 
ficient amplitude to allow the formation of low mass stars 
or even brown dwarfs within complex layers of post-shock 
gas, the turbulence may directly affect the mass distribu- 
tion of pre-stellar cores and stars (p^adoan & N ordlundl l2002[ 
iPadoan et albOOTtlHennebeUe & Chabrieil20 08. 2009). Fur- 
thermore, the turbulence must be at least partly responsible 
for the low star formation rate per free-fall time observed in 
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most environments (iKrumholz & Tanll2007h . because the tur- 
bulent energy generally exceeds the gravitational energy on 
small scales within molecular clouds (the vi rial parameter is 
almost alw ays larger than unity, as sh own by iFalgarone et al.l 
(Il992h and iRosolowsky et alJ (l2008h in Perseus). Theoreti- 
cal models of the star formation rate b ased on the effect of 
turbu l ence have recently been pro posed (IKrumholz & McKe3 
l200llPadoan & Nordlundl201lb . 

The importance of turbulence in the process of star for- 
mation provides an opportunity for theoretical modeling, be- 
cause one can assume that molecular clouds follow the uni- 
versal statistics of turbulent flows, for example with respect to 
the probability density function (PDF) of gas density and the 
scaling of velocity differences. The turbulence is also a chal- 
lenge for numerical simulations of star formation, because the 
limited dynamical range of the simulations cannot always ap- 
proximate well enough the scale-free behavior of the turbu- 
lent flow. The Kolmogorov dissipation scale, is the small- 
est turbulent scale below which viscous dissipation becomes 
dominant. It can be computed as 7]k = (^^''/e)'''"*, where v is 
the kinematic viscosity and e the mean dissipation rate of the 
turbulence. The kinematic viscosity can be approximated as 

V ~ Vfnjian), where i^th is the gas thermal velocity, n is the 
gas mean number density, and cr « 5 x 10"'^ cm- is the gas 
collisional cross section. The mean dissipation rate can be 
estimated as e ^ ji, where ^ is a scale within the inertial 
range of the turbulence, and v is the rms velocity at the scale 
i. In molecular clouds, assuming the lLarsonl(II981l) relations 

V ^ \ km s"'(i?/lpc)""*- and n ^ lO-' cm"-'(£/Ipc)~', a gas 
temperature of 10 K, and a driving scale of ~ 70 pc, we obtain 
77k lO''* cm, well below the characteristic spatial resolution 
of the gas dynamics in star formation simulations. 

The dynamic range limitation of the simulations can be ex- 
pressed in terms of the Reynolds number The Reynolds num- 
ber estimates the relative importance of the nonlinear advec- 
tion term and the viscosity term in the Navier-Stokes equa- 
tion. Re = u„nsi3/i/, where Uims = \/ (f^) is the flow rms ve- 
locity, a is the integral scale of the turbulence (of the order 
of the energy injection scale). The Reynolds number can also 
be expressed as Re = {Lji'iY^'^l^ . Based on the same assump- 
tions used above to derive ?7k, we obtain Re 10^ for typical 
molecular cloud values. At present, the largest simulations 
of supersonic turb ulence may achieve an effective Reyn olds 
number Re - 10"* (iKritsuk et alJl2009aLlJones et alJl201Ih . 

Numerical simulations are incapable of describing the 
smallest structures of magnetic fields in star-forming clouds. 
The characteristic magnetic diffusivity, r/, of the cold inter- 
stellar gas is much smaller than the kinematic viscosity, v. As 
a result, magnetic fields can develop complex structures on 
scales much smaller than the Kolmogorov dissipation scale, 
77k, where the velocity field is smooth. Introducing the mag- 
netic Prandtl number, Pm, defined as the ratio of viscosity and 
diffusivity, Pm = v/rj, this regime is characterized by the con- 
dition Pm ^ 1. The magnetic diffusivity can be expressed as 
■q = c^mei'en/47rnee^ (cgs), where c is the speed of light, the 
electron mass, the collision frequency of electrons with 
neutrals, «e the number density of electrons, and e the elec- 
tron charge. This expression neglects electron-ion collisions, 
because at the low ionization fractions and temperatures of 
molecular clouds the dominant friction force on the electrons 
is from collisions with neutrals. The collision frequency of 
electrons with neutrals can be written as i/gn = «nCfth,e, where 
ria is the number density of neutrals (~ n in molecular clouds). 



(T is the gas collision cross section given above, and Uth.e is the 
thermal velocity of the electrons. The magnetic Prandtl num- 
ber is then given by Pm w 2 x 10^(xi/10"^)(/?/1000cm-^)-', 
where Xi is the ionization fraction. 

Numerical simulations without explicit viscosity and mag- 
netic diffusivity usually have effective values of Pm ^ 1, very 
far from the conditions in molecular clouds. If the magnetic 
field strength is determined self-consistently by a small-scale 
turbulent dynamo, this numerical limitation may cause an ar- 
tificially low magnetic field strength in low-resolution simula- 
tions, or in simulations based on MHD solvers with large ef- 
fective magnetic diffusivity. Such simulations may not reach 
the critical value of the magnetic Reynolds number, Rm, re- 
quired by the turbulent dynamo. The magnetic Reynolds 
number is defined as Rm = RePm = Wi-ms'C / rj. Its critical value 
for the turbulent dynamo in supersonic turbulence was found 
to be Rmci-it ~ 80 in the regime with Pm ~ 1 and for a sonic 
rms Mach number Ms ~ 2.5, where Ms = Vymsl cs is the ratio 
of the flow rms velocity and th e speed of sound (iHaugen et al.l 
12004 . iFederrath et"aLl (l20IIh find Rmcrit « 40 for transonic 
turbulence, driven by the gravitational collapse of a dense, 
magnetized gas cloud. 

Besides the effective Re and Pm, the other two non- 
dimensional parameters of isothermal MHD turbulent simu- 
lations are the rms sonic Mach number, defined above, and 
the rms Alfvenic Mach number, Mq.a = Urms/uo,A, where uq.a 
is the mean Alfven speed defined as wq.a = Bq/ yJAiipQ, and 
Bfi and po are the mean magnetic field and mean gas density, 
respectively. The initial conditions of the numerical test de- 
scribed in this work have Ms ~ 9 and Mq.a ~ 30. In the test 
runs, the value of Ms decreases with time (no driving force is 
used), as shown in the left panels of Figure [Tj vq^a is instead 
constant, because both Z?o and po are conserved quantities in 
the simulations. However, the rms value of the magnetic field 
strength, ^/(W), depends on both Bq and Wims . In these sim- 
ulations Bo is very low, and the turbulence is highly super- 
Alfvenic, meaning that Wi-„is ^ I'o.A- In this regime, the mag- 
netic field is locally amplified by compression and stretching 
resulting in a statistically steady state with ^ {B~) ^ Bq. The 
rms Alfvenic Mach number defined in terms of the mean mag- 
netic and kinetic energies. Ma = \J {pv'^) / (B^/47r) « 4.4 and 
decreases with time as the turbulence decays, as shown by the 
right panels of Figure [T] 

Based on the observed dependence of the v elocity disper- 
sion on spatial scale in molecular clouds (e.g.. lLarsonllI98lt 
iHeyer & Bruntll2004l) . the initial value of Ms in our test runs 
is relevant to star-forming regions on scales of a few par- 
sees. The super- Alfve nic nature of molecu l ar clo ud turbu- 
lence was suggested by iPadoan & Nordlundl ( II999I). and has 
received further support in more recent work (iLunttilaet al.l 
l2008ll2009l:IPadoan et al.llMol: IKritsuk etani20I Ih . 

One way to assess the ability of numerical simulations to 
approximate the behavior of turbulent flows is to study the 
power spectra of relevant quantities, such as velocity and 
magnetic fields. The interpretation of velocity power spec- 
tra from numerical simulations face the following problems: 
(1) the limited extent of the inertial range of turbulence due 
to the limited range of spatial scales discussed above (or even 
the complete absence of an inertial range in the case of low 
resolution simulations); (2) the emergence of the bottleneck 
effect in hydrodynamic simulations (e.g.. jFalkovichl 1 1 9941: 
iDobler et al.ll2003l: IHaugen & Brandenburdl2004l) as soon as 
the numerical resolution is large enough to generate an iner- 
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tial range; (3) the dependence of the power spectrum on the 
numerical schemes; (4) the dependence of the numerical res- 
olution necessary for convergence on the numerical method. 

This work addresses the above problems and the general 
issue of the quality of MHD codes with respect to the de- 
scription of highly supersonic and super- Alfvenic isothermal 
turbulent flows. We do not study the quality of simulations 
of the gravitational collapse of gravitationally unstable re- 
gions with adaptive mesh refinement or Lagrangian meth- 
ods in this paper Although most star formation simulations 
eventually take advantage of such techniques, here we fo- 
cus on the simulations of turbulent flows where gravity is 
neglected. This work considers high-resolution simulations 
of MHD turbulence, while related stu dies of nonmagnetized 
flo ws have b een recently published bv lKitsionas et al.l ( l2009h 
and lPrice & Fe derrath (20Tol). 

The paper is organized as follows. In Section 2, we de- 
scribe the simulation setup. In Section 3, we introduce the 
algorithms used. In Section 4, we discuss the diagnostic tech- 
niques utilized in the paper In Section 5 we present the re- 
sults from each code, and in Section 6 we discuss the impact 
of method design on the numerical dissipation properties. Fi- 
nally, Section 7 summarizes our conclusions. 

2. THE TURBULENCE DECAY TEST PROBLEM 

Modern numerical methods for astrophysical turbulence 
simulations are designed to produce approximations to the 
limit of viscous and resistive solutions as the viscosity and 
magnetic diffusivity are reduced to zero. Numerical experi- 
ments carried out with such methods can be viewed as im - 

t'icit large eddy simulations, or ILES (iGrinstein et al.ll2007h . 
/tine et al. ( 2000) de monst rated that Euler solvers, like PPM 
:olella & Woodwardl[T984h . are more efficient than Navier- 
Stokes solvers in providi ng a better scale se paration at a given 
grid resolution (see also lBenzi et al.ll2b08l) . Here we employ 
the same ILES technique for MHD simulations of decaying 
supersonic turbulence. The numerical methods we compare 
differ in their implicit subgrid models and the focus of this pa- 
per is on understanding the origin of those differences, which 
could help to improve our methods. 

We, thus, solve numerically the system of MHD equations 
for an ideal isothermal gas in a cubic domain of size L with 
periodic boundary conditions: 



dp 
dt 



dpu 
~di' 



+ V- 



+ V-(pu) = 0, 



PUU-BB+ p + — I 



dB 

— + V (uB-Bu) = 0. 

at 



(1) 



(2) 



(3) 



Here, p and u are the gas density and velocity, B is the mag- 
netic field strength, p is the gas pressure, and I is the unit 
tensor 

All numerical methods discussed in this paper are designed 
to conserve mass, momentum and magnetic flux, and attempt 
to keep V • B = to the machine precision. All methods are 
formulated to approximate the ideal MHD Equations (l)-(3). 
However, due to the finite numerical viscosity and magnetic 
diffusivity, as well as artificial viscosity and diffusivity added 
for numerical concerns, the actual equations evolved will have 
additional dissipation terms on (2) and (3). The exact nature 
of these dissipation terms is method-dependent. 



In this section and below, we use dimensionless code units, 
such that the domain size L = 1 ; the gas density p is given in 
units of the mean gas density po; the gas pressure p is given 
in the units of the uniform initial pressure po, and the velocity 
u is given in units of the sound speed, u = v/c^. The uniform 
mean magnetic field is Bq = yj2/ = 0.3, where the ratio of 
thermal-to-magnetic pressure /3o = 22. The code units also 
imply that B incorporates the 1 /Att factor so that the magnetic 
pressure is given by B^/2 in the code units. 

Initial conditions for the decay test were generated in 2007 
with an earlier, non-conservative version of the STAGGER 
code on a lOOO^* grid using a time-dependent random large- 
scale (k/kuim < 2, where k^m = I-k/L) isotropic solenoidal 
force (acceleration) F to stir the gas and reach an rms sonic 
Mach number Msf) ~ 9. There was no forcing in the induc- 
tion Equation (3), so the rms magnetic field was passively 
amplified through interaction with the velocity field. The 
model was initiated with a uniform density po and pressure 
Po, random large-scale velocity field Uo, and a uniform mag- 
netic field Bo aligned with the z-coordinate direction. To 
achieve a saturated turbulent state, the flow was evolved with 
the STAGGER code for three dynamical times (defined as 
fd = L/2Ms,o). Assuming an initial Ms.o = 10, fd = 0.05 in the 
code units determined by the box sound crossing time. In the 
saturated turbulent state, the level of magnetic fluctuations is 
~ 50 times higher than Bq, i.e., B = Bo + b, where /jnns ^ 
and (b) = 0. 

The actual test runs were performed at grid resolutions of 
256^ 512^ and in a few cases (PPML and ZEUS) also 1024^ 
cells. Data regridding utilized conservative interpolation of 
hydrodynamic variables while a vanishing V • B in the inter- 
polated initial states was enforced with V • B cleaning. The 
evolution of decaying turbulence (F = 0) was followed for 
Af = 0.2 = 4fd and 10 flow snapshots equally spaced in time 
were recorded for subsequent analysis. The timing of these 
snapshots in the adopted time units is as follows: ti = 0.02, 
t2 = 0.04, . . . , fio = 0.2, assuming f = corresponds to the end 
of the initial forcing period. 

3. NUMERICAL METHODS AND IMPLEMENTATIONS 
3.L ENZO2.0 



ENZO's dO'Sheaet al.ll2005l) MHD scheme dWang & Abell 
|2009|) employs the following components: second- order spa- 
tial in terpolation via the Piecewise Linear Method (Ivan Leed 
Il979l) : second-order time integration v ia a second-order 
Runge-Kutta method dShu & Osheill988 '); the HLL Riemann 
solver f or computation of int erface fluxes (Harten et al. 1983); 
and the iDedner et al.l (l2002l) scheme for maintaining the di- 
vergence of the magnetic field close to zero. The code is for- 
mally second-order accurate in both time and space. These 
one-dimensional components are combined to form a three- 
dimensional method in a directionally unsplit manner, with 
the Runge-Kutta integration mediating the wave information 
between the three flux computations. The slope limiter 9, 
which contr ols the sharpness of the reconstruction, was set 
at 1.5 as in dWang et al.ll2008h . Larger values were tried for 
256^ grids, without significant change to the solution. 

3.2. FLASH 3 

The FLASH3 dFrvxell et al.ll2000t iDubey et alJlIOOSh sim- 
ulations presented in this study hav e used a completel y 
new MHD scheme implementation dLee & Deand l2009h . 
The solver adopts a dimensionally unsplit integration on 
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a staggered grid (Unsplit Staggered Mesh), for the muhi- 
dimensional MHD formulation, based on a finite-volume, 
higher-order Godunov method. A new second-order data 
reconstruction-evolution method, ex tended f rom t he corner 
transport upwind (CTU) approach of lColellal (1 1990b has been 
used, which guarantees proper evolution of in-plane dynamics 
of magnetic fields. The importance of the in-plane field evolu- 
tion is described and test ed in the field-loop advection test in 
iGardiner & Stonel (120051) . The unsplit staggered mesh solver 
(USM) has also shown a successful performan ce on this test, 
maintaining a correct in-plane field dynamics (iLee & Deanel 
120091) . The algorithm uses a new "multidimensional charac- 
teristics analysis" in calculating transverse fluxes. This ap- 
proach is advantageous and efficient because it does not re- 
quire solving a set of Riemann problems for updating trans- 
verse fluxes. High Mach number turbulent flows require a 
precise and positive-preserving solver capable of resolving 
complex shock structures while keeping numerical diffusion 
as sma U as possible. We therefor e chose the HLLD Riemann 
solver (iMiyoshi & Kusanoll2005h . which greatly improves the 
robustness and accuracy of supersonic MHD turbulence simu- 
lations as the Roe solver easily fails to preserve positive states 
of density and/or pressure in strong rarefaction waves. For 
further enhancing solution accuracy and stability, we chose 
a hybrid limiter that uses the compressive van Leer's slope 
limiter for linearly degenerate waves and the more diffusive 
minmod limiter for genuinely nonlinear waves. 

3.3. KT-MHD 

The KT-MHD code is an implementation of a semidiscrete 
centra l-difference scheme developed by Kurganov & Tadmo^ 
(I2OOOI) . The total time derivative of the hydrodynamic quanti- 
ties is computed using the flux definition of the Kurganov- 
Tadmor scheme, a higher order extension of the Lax- 
Friedrichs scheme. The flux values are evaluated at the cell 
interfaces. The corresponding point values of the conserved 
quantities are interpolated to the cell interfaces via a third- 
ord er CWENO scheme in thr ee space dimensions follow- 
ing IB albas and Tadmo^ (l2006h . The averages of the mag- 
netic field components reside at the cell interfaces and are 
reconstructed in diagonal direction, also using a third-order 
CWENO scheme. The smoothness indicators (and thereby 
the nonlinear weights) of the CWENO scheme are based on 
the density field, only. Componentswise smoothness indica- 
tors have shown to lead to a much higher numerical viscosity. 
The total time derivative of the magnetic fi eld is compute d 
by a constrained transport (CT) method of IZieglerl (l2004b . 
The resulting set of ordinary differential equations is inte- 
grated in time by a fourth-order Runge-Kutta scheme. The 
code uses a regular grid and the so-called pencil decompo- 
sition in its MPI-parallel implementation. The idea of com- 
bining the Kurganov-Tadmor central-difference scheme with 
a CT method for the magnetic field update was first imple- 
mented by Ralf Kissmann and published in his PhD the sis 
(lKissmannll2006t) and is used bv lDreher and Grauerl (l2005l) in 
their Racoon code. 

3.4. LL-MHD 

The CT-based LL-MHD solver (ICoUins etalj|2010l) em- 
pl oys a div e rgence preserving higher-order Godunov method 
of iLi et all (l2008h . which uses second-order spatial recon- 
struction and second order time reconstruction to compute the 
interface states, and the isothermal HLLD Riemann solver of 
iMignonj ( l2007h to compute the flux from those reconstructed 



states. This is done in a directionally split fashion, with the 
order permutation of Strang to preserve th e second-order ac- 
curacy . The solver uses the CT method of IGardiner & Stonel 
( I2OO5I) to maintain the divergence-free evolution of the mag- 
netic field. LL-MHD is also installed in the AMR code 
ENZO, and has been used to study a range of astr ophysi- 
cal phenome na, from galaxy clu sters (IXu et al.ll2010h to pre- 
stellar cores (IColUns et al.ll20TTh . 

3.5. PLUTO 3.1 

The PLUTO code (iMignond l2009l) is a highly modular, 
multi-dimensional and multi-geometry code that can be ap- 
plied to relativistic or non-relativistic MHD or HD (hydro- 
dynamic) flows. PLUTO comprises several numerical meth- 
ods, like the high-order conserva tive finite-difference diver- 
gence cleaning MHD method CMignone & Tzeferacosll201C)l) 
as well as finite-volume CTU schemes jMignone et al.B2010l) . 
The latest version of the PLUTO code (V. 3.1—2010 Au- 
gust) allows to choose between several space reconstruction 
and time integration methods as well as several approximate 
Riemann solvers including HLL, HLLC, HLLD or the Roe 
Riemann solver For the MHD forr nulation one can ch oose 
between the eight-wave for mulation jPowell et al.ll 19991) . the 
divergence cleaning method ( iDedner et al.ll2002h . and the CT 
method. The possibility to switch between several numer- 
ical methods allows to handle a wide range of astrophysi- 
cal problems. For this test we used the accurate Roe Rie- 
mann solver in corabinat ion with a third-order reconstruction 
(ICada & Torrilhonll2009l) . characteristic variable limiting, the 
Runge-Kutta 3 time integration and the iPowell et alj d 19991) 
eight-wave MHD formalism; three-dimensional effects were 
incorporated by way of the Runge-Kutta integration, without 
the use of the transverse flux gradients used in CTU. Courant 
number was set to 0.3. 

3.6. PPML 

The piecewise para bolic method on a local stencil (PPML, 
lUstyugov et al.ll2009l) is a compact stencil variant o f the pop- 
ular PPM algorithm (IColeUa & WoodwardI 11984 for com- 
pressible magnetohydrodynamics. The principal difference 
between PPML and PPM is that cell interface states are 
evolved rather that reconstructed at every time step, result- 
ing in a more compact stencil. The interface states are 
evolved using Riemann invariants containing all transverse 
derivative information. The conservation laws are updated 
in an unsplit fashion, making the scheme fully multidimen- 
sional. Divergence-free evolution of the magnetic field is 
maintained using t he hig her-order-accurate CT technique of 
IGardiner & Stonel (l2005l) . The method employs monotonic- 
ity constraints to preserve the order of scheme in points of 
local extrem a (Suresh and Huynh 1997; Balsara & Shu 2000t 
iRider et aHIToOTb . To preser v e mon otonicity in multidimen- 
sions a method from iBarthI (1 19901) is additionally applied. 
An updated component of the electric field at a cell bound- 
ary is calculated by averaging the quantities obtained from 
known compone nts of flux-vectors and va lues of gradient of 
the electric field (IGardiner & Stonell2005l) . The performance 
of PPML was tested on several numerical problems, which 
demonstrated it s high accuracy on bot h smooth and discontin- 
uous solutions dUstvugov et al.ll2"009l) . Simulations of super- 
sonic magnetized turbulence in three dimensions with PPML 
show that low dissipation and wide spectral bandwidth of this 
method make it an ideal can didate for direct turbulence simu- 
lations ( iKritsuk et al.ll2009allbh . 
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3.7. RAMSES 

RAMSES jTevssied l2002h is an unsplit Godunov AMR 
scheme with a second-order total variation diminishing spa- 
tial reconstruction using the Monotonized Central slope lim- 
iter Magnetic field is updated using the CT method, using 
two-dimensional Riemann problem at cell edges to compute 
the electro-motive-force that enters into the induction equa- 
tion. The magnetic field divergence, expressed in an inte- 
gral form on cell faces, is therefore zero down to machine 
accuracy. Conservative variables are updated by solving one- 
dimensional Riemann problems at cell faces. Both the one- 
dimensional and the two-dimensional Riemann solvers are 
based on the H LLD MHD approximate Riemann solution 
dMivoshi & Kusano 2005) . More details on the MHD scheme 
can be found in Teyssier et al. (2006) and Fromang et al. 
(2006). 

3.8. STAGGER 

The STAGGER Code is originally based on a code 
develop ed as part of the PhD thesis of Klaus Gals- 
gaard (iGalsgaardI Il996l) . Several versions exist, and 
the code is used i n many different circumstances 



1998. 2000; ,Stein & Nordlundl 


19981 lAsDlund et all 


20001 iPadoan et all 120041: iGudiksen & Nordlundl I2005t 


Braithwaite & Nordlund 2006 


; Archontis et al. l2007l 


Lunttila et al.l 120091 IStein et al. 


1201 ll IPadoan & Nordlundl 


201 1|). 



In the context of solar and stellar physics it is equipped with 
a multi-frequency radiative transfer module and a comprehen- 
sive equation of state module that includes a large number 
of atomic and molecular species, to be able to compute re- 
alistic three-dimensional models of the near-surface layers of 
stars. The widths, shifts, and asymmetries of synthetic spec- 
tral lines computed from such models exemplifies some of the 
most precise agreements between three-dimen sional numeri- 
cal sim ulations and astrophysical observations (lAsplund et alj 
l2000h . 

In t he context of supersonic turbu lence studies earlier 
works (IPadoan et al.ll 1 997l 1 1 998l l2000l) were based on a non- 
conservative version of the code, which evolved the primi- 
tive variables Inp, u, and B. The "per-unit-mass" formulation 
based on these variables is simple and robust, but has the dis- 
advantage that mass and momentum are not conserved exactly 
by the discretized equations. 

The current version of the code instead uses the per-volume 
variables p, pu, and pE, where E is the internal energy per unit 
mass, allowing a discretization that explicitly conserves mass, 
momentum, energy, and magnetic flux. In the isothermal case 
of relevance here the code solves these equations: 



dt 
dpu 

'dt 



= -V -(puu + r)- Vp+(V X B) X B, 
= -V X (-U xB + TyV xB), 



where t is the viscous stress tensor, which we write as 

Tij=-pvSij, 

and Sij is the strain rate 

^ 1 / duj ^ du. 



(4) 



(5) 



(6) 



(7) 



(8) 



The viscosity v and magnetic diffusivity ?y are spatially de- 
pendent, in a manner reminiscent of the Richtmyer & Morton 
formulation, with 

= («lMW + «2'^M^)Ai , (9) 

where Ai is the mesh size, «w is the wave speed, 5u^ is the 
positive part of a second order approximation of -As V • u. 
The magnetic diffusivity is taken to be 

ri = nB{n\Ui^ + n25ug)As , (10) 

where 5u^ is analogous to 5u^, except only the component of 
the velocity perpendicular to B is counted. 

Here, ni, 112, and are numerical coefficients of the order 
of unity. The «iMw term, where «i ^ 0.03 is a relatively small 
constant, is needed to provide stabilization and a weak dis- 
persion of linear waves, while the n25u" term, with «2 ^ 0.5, 
provides enhanced dissipation in shocks, where the rate of 
convergence -V • u is large. Since the magnetic field is in- 
sensitive to motions parallel to the field, only perpendicular 
motions are gauged by the corresponding magnetic diffusiv- 
ity term, hb is essentially an inverse magnetic Prandtl number 

The tensor formulation of the viscosity ensures that the vis- 
cous force is insensitive to the coordinate system orientation, 
thereby avoiding artificial grid-alignment. 

3.9. ZEUS-MP 

ZEUS-MP is a widely used, multiphysics, massively par- 
allel, message-passing Eulerian code for astrophysical fluid 
dynamic simulations in three dimensions. ZEUS-MP is a dis- 
tributed memory version of the shared-memory code ZEUS- 
3D that uses block domain decomposition to achieve scal- 
able parallelism. The code includes hydrodynamics, mag- 
netohydrodynamics, and self-gravity. The HD and MHD al- 
gorithms are bas ed on the method of fini te differences on a 
staggered mesh (IStone & Normanlll992allbl) . which incorpo- 
rates a second-o rder-accurate, monotonic advection scheme 
dvan Leeilll977h . The MHD algorithm is suited for multidi- 
mensional flows using the me thod of character i stics s cheme 
(MOC-CT) first suggested by lHawley & Stond (Il995h . Ad- 
vection is performed in a series of directional sweeps that 
are cyclically permuted at each time step. Because ZEUS- 
MP is designed for large simulations on parallel comput- 
ing platforms, considerable attention is paid to the paral- 
lel performance characteristics of each module in the code. 
Complet e discussion on all algorithms in ZEUS-MP can be 
found in lHayes et"an (l2006l) . All the MHD turbulence decay 
simulations performed using ZEUS-MP in this paper use a 
quadratic (von Neumann-Richtmyer) artificial viscosity coef- 
ficient qcon of 2.0 and a Courant number of 0.5. 

4. DATA ANALYSIS 
4. 1 . Power Spectra 

Given a vector field u(r) discretized on a mesh i,j,k with 
Ui,j,k one can compute a power spectrum from the three- 
dimensional Fourier transform u,. by summing the magni- 
tudes squared, |u/.jjt|-, overk-shells with K„ < |k, < K„+i. 
If the Fourier transform coefficients Uij.k are normalized so 
the rms value of the corresponding function in real space is 
equal to unity, then the sum of the squares in Fourier space 
is equal to the average of the function squared in real space 
(Parseval's relation): 



rms 



i.j,k 



"m>I' = ]5^EI"'.mI', (11) 
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TABLE 1 

Selected Numeric Values for the Decay Test 









2rj + 4/3A'= " 


u-bandwidth B-bandwidth ' 


otir 1 not . ~\ s 


ENZO 


LOOl 


0.78 


0.93 


0.92 


0.19 


0.07 


0.60 


FLASH 


LOGO 


0.94 


0.85 


1.38 


0.15 


0.20 


0.27 


KT-MHD 


L041 


0.85 


0.89 


1.30 


0.20 


0.13 


0.86 


LL-MHD 


L062 


0.81 


1.02 


0.80 


0.22 


0.10 


0.29 


PLUTO 


L077 


0.92 


1.03 


1.14 


0.20 


0.12 


0.32 


PPML 


L043 


0.92 


1.20 


1.46 


0.24 


0.20 


0.32 


RAMSES 


L069 


0.87 


1.07 


1.18 


0.24 


0.09 


0.33 


STAGGER 


LOOS 


0.70 


1.93 


0.79 


0.28 


0.07 


0.31 


ZEUS 


L037 


0.83 


0.76 


1.01 


0.16 


0.10 


0.27 



" Mean specific kinetic energy density at / = 0.2 normalized by tlie reference solution; see Section 5.1 and Figure 1. 
^ Mean magnetic energy density alt = 0.2 nomialized by the reference solution; see Section 5.1 and Figure 1. 

A proxy for the mean dissipation rate of specific kinetic energy at f = 0.02; see Section 5.2 and Figure 2 left. 

A proxy for the mean dissipation rate of magnetic energy axt = 0.02; see Section 5.2 and Figure 2 right. 

Effective spectral bandwidth for the velocity; see Section 5.4 and Figure 4 left. 
' Effective spectral bandwidth for the magnetic field; see Section 5.4 and Figure 4 right. 

i Ratio of dilatational-to-solenoidal power averaged over k/k„^\„ > 100 at ^ = 0.2; see Section 5.5 and Figure 5 right. 



where is the total number of /, j,k points. 

In the codes used to analyze the results for the current pa- 
per we use the real valued Fast Fourier Transform routine 
sr f f t f from the f f tpack software package, which returns 
coefficients af:=N/2 for sine and cosine functions, except for 
the DC and Nyquist components, which return coefficients 
A^. Proper power normalization requires that sine and cosine 
components contribute power 1/2, and the returned coeffi- 
cients should thus be multiplied with Vl/N, except for the 
DC and Nyquist components (which are the first and last coef- 
ficients returned from srf f tf), which should be multiplied 
with 1 /N. 

The power spectrum P{k) expresses how much of the power 
falls in each fc-interval. If the power is collected in discrete 
bins, 

Pn= l"'.^-.'^!'' (12) 

K„<\kiji\<K„+i 

then the total power can also be expressed as 

rms- = ^P,„ (13) 

n 

where the sum is taken over all bins. 

To illustrate the power spectrum P„ graphically one needs 
to assign a wavenumber k„ to each bin. A natural but not 
quite optimal choice is to use the midpoint of the bin; k„ = 
(K„+K„+\)/2. A better choice is to use the mean of the 
wavenumbers that actually fall inside the bin (to see why this 
is better consider a case with very wide bins and a function 
with power at only a few discrete wavenumbers). 

It turns out that one gets smoother power spectra if one as- 
signs a value 

P'^J-^iKl,-K^,)^Y.\^,j,\' (14) 
bi" bin 

rather than 

^« = $]|u,V,-i|' (15) 
bin 

to each bin. In other words: power spectra (at least those mea- 
suring fluid flow properties) become smoother if they measure 
the average power in a shell (times the shell volume) rather 



than the total power One can interpret this to mean that fluid 
flow properties are encoded in Fourier amplitudes as a func- 
tion of wavenumber, rather than in total power of Fourier am- 
plitudes in a shell. If (and only if) this is the case then the 
power spectrum fluctuates (as observed) down (or up) if by 
chance a shell contains fewer (or more) discrete wavenum- 
bers than expected. 

To be able to recover P„ from P,' (e.g. for use in Parse- 
val's relation) it is necessary to record the number of discrete 
Fourier amplitudes in each bin, N-^\^ in Equation (fl4l i above. 

Note also, that in order for Parseval's relation to be exact 
for three-dimensional power spectra, all Fourier components 
need to be included, which means that the A:-scale should re- 
ally extend to a maximum value of -v/S^n, where is the 
one-dimensional Nyquist frequency. Nevertheless, here we 
follow the common practice to truncate the three-dimensional 
power spectra at k^^. 

4.2. Helmholtz Projections 
The divergence of a vector field f, j t, with Fourier transform 
coefficients tij.k, is 

FT(V-f) = ik,v,A.-f, (16) 

The vector coefficients f, y ^ may be split into a component 
parallel to k, and a remaining component, which is perpen- 
dicular to k/.jjt: 

fl^, = ^iJ,,{^^iJ,t\j,k)/\^^iJA^ (i7) 

and 

^j± = l.j^-^jr (18) 

Taking the divergence of the latter, we have 

i%j,k~iljk = i%.j.k ■h.j.k-i%.j.k ■'^i.j.ki^.j.k ■~^i.j.k) /%.j.k? = 0. 

(19) 

The inverse transform based on the tj-j ^ coefficients is thus 

solenoidal, while the inverse transform based on f|'^- is purely 
compressional. 
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Fig. 1. — Time evolution of tlie mean specific kinetic energy (top left), magnetic energy (top right), as well as sonic (bottom left), and Alfvenic (bottom 
right) rms Mach numbers at grid resolution of 512' cells. Note, that the kinetic energy and sonic Mach number are rather insensitive to the details of numerical 
dissipation, while the evolution of magnetic energy and Alfvenic Mach number display significant dependence on the numerical magnetic diffusivity. 




Fig. 2. — Time evolution of — eKRseff = 2f2 4-4/3A defined in Section \5?2\ (left) and the mean-squared current, 7^ = (| V X B|^) (right). These "small-scale" 
measures of turbulent fluctuations are sensitive to the details of numerical diffusivity and highhght differences between the methods. 



5. RESULTS 

There is a great variety of interesting statistical measures in 
magnetized supersonic turbulent flows to study and compare. 
The KITP07 project originally implied a comparison of den- 
sity structures in physical space using projections and slices, 
PDFs of the gas density, various power spectra and structure 



functions, time evolution of some global and local average 
quantities, etc. Most of these data provided by individual con- 
tributors to the project can be accessed electronically via the 
wiki Web site at KITP under the rubric Star Formation Test 
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Problems. 

In this paper we mainly focus on the statistics of the ba- 
sic MHD fields (the so-called primitive variables) since those 
are easier to interpret and link back to the essential features 
of the numerical methods. Since the system of equations and 
the initial and boundary conditions are the same for all codes, 
the only source of differences in the numerical solutions is 
numerical dissipation. In this section, we discuss the sensi- 
tivity of various turbulence diagnostics to the numerics and 
describe a set of statistical measures that allow us to assess 
the quality of different algorithms we compare. We begin 
with global averages over the periodic domain and then con- 
tinue with analysis of power spectra. We avoid discussion of 
density statistics, even though these are important for numer- 
ical star formation studies. That discussion would be more 
appropriate in a context of driven turbulence, where time- 
averages over many flow snapshots help to red uce the strong 
statistical no ise associated with the density (e.g.. lKritsuk et alj 
I2OO6I l2007h . The density statistics in a similar context have 
been discussed in detail elsewhere (e.g.. lKitsionas et ani2009l : 
Iftice & Federrath..20T0l) . 

5.1. Mean Kinetic and Magnetic Energy, mis Mach Numbers 

The evolution of the mean kinetic energy is captured per- 
fectly well by all methods, except for some small (< 8% by 
f = 0.2, see Table 1) differences that become noticeable at 
t > 0.12 in Figure [T] This particular quantity is known to 
converge rather early in nonmagnetized compressible turbu- 
lence simulations, and the same is true for super- Alfvenic tur- 
bulence (Lemaster& Stone 2009; Kiitsuket al. 2009b). The 
velocity power spectrum P{u,k) ^ k" has an inertial range 
slope a e [-2,-5/3] that depends on the sonic Mach num- 
ber (see, for instance, Figure |3] below). Because the spec- 
tral slope is so steep, the mean specific kinetic energy den- 
sity, EjQ = (u~/2) = P{u,k)dk/2, is strongly dominated by 
large scales. If the resolution is sufficient to properly capture 
the structure of large scale flow in the computational domain, 
the energy convergence is achieved. This is apparently the 
case in our 512^ simulations, see the left panels of Figure [T] 
For the 256^ model (not shown), the result is very similar. 
We thus conclude that in the decaying turbulence problem the 
mean kinetic energy is not sensitive to variations in small- 
scale numerical diffusivity between different methods. 

The mean magnetic energy density, Em = (B^/2) = 

P{B,k)dk/2, appears to be more sensitive to variations 
in small-scale numerical kinetic and magnetic diffusivity in 
super- Alfvenic simulations, see Figure [T] right panels, and 
Table 1. Most of the methods show an early-time increase 
in magnetic energy, but asymptotically, after saturation is 
reached, all of them show very similar decay rates Em/Em- 
The saturated level of Em in the initial flow snapshot gener- 
ated with the original, non-conservative version of the STAG- 
GER code is lower than most other codes would produce, 
except for perhaps LL-MHD and ENZO. To compensate for 
this deficiency of magnetic energy in the initial conditions, 
FLASH and PPML add about 7% to the initial Em by the time 
fi = 0.02, when the first flow snapshot is recorded. KT-MHD 
increases Em by - 4%, PLUTO, ZEUS and RAMSES add 
- l%-2%. The level of reached with the old STAGGER 
code is roughly consistent with that of ENZO and LL-MHD. 
The new, conservative STAGGER appears to be more diffu- 

http://kitpstarfoiTnation07.wikispaces.com/Star+Formation+Test+Problems 



sive than all other methods as far as the magnetic energy den- 
sity is concerned. 

To understand why the magnetic energy levels are differ- 
ent, one needs to recall that the convergence rate for Em with 
the grid resolution is rather slow at 512''. For instance, with 
PPML, the saturated levels of Em in driven simulations con- 
tinuously grow as grid resolution improves an d the conver- 
gence is ex pected only at 2048^* or even higher (iKritsuk et al.l 
2009b; Jon es et al.l 120111) .'^ The slow convergence of Em in 
super-Alfvenic runs is not surprising because of their rather 
shallow magnetic energy spectra (see Figure |3] below and 
note that the magnetic spectra are plotted noncompensated). 
In such weakly magnetized (Z?o <C ^ims) isotropic nonhelical 
flows, turbulence amplifies the rms magnetic field fluctuations 
by stretching and tangling the field lines primarily on small 
scales until an equilibrium is reached between the rms field 
amplification and dissipation. The saturated level of Em for 
a given sonic Mach number, Ms, and mean magnetic field 
strength, Bq, would naturally depend on the effective mag- 
netic Prandtl number and on the effective magnetic Reynolds 
number We discuss the relative standing of the methods in 
terms of Rm and Pm in the next section and show that the sat- 
uration level of magnetic energy indeed correlates with Rm 
and Pm. Overall, by f = 0.2, the deviations of from the 
reference solution defined in Section 5.4 can be as large as 
30%, see Table 1. 

Turbulence regimes simulated with different methods dif- 
fer slightly in their rms Alfvenic Mach numbers. Ma = 
■\/2 (pu^) I (-6^). We use this proxy for the Alfvenic Mach 
number instead of ^J2{pu-/B~) since in the latter case the lo- 
cations where B (nearly) vanishes would produce arbitrarily 
large contributions to the mean making this measure unstable. 

At f = 0.1, the least super-Alfvenic regime is achieved with 
FLASH and PPML, followed by PLUTO, KT-MHD, RAM- 
SES, ZEUS, ENZO, LL-MHD, and the new STAGGER in the 
order of increasing Ma. Note, that the ranking of the methods 
is essentially the same as for Em, which, unlike Ma, does not 
depend on the gas density. This similarity can be explained 
by a limited sensitivity of our chosen proxy for Ma to corre- 
lations betwe en density and veloci ty or between density and 
field strength (iKritsuk et al.ll2009ah . 

Finally note that while Zsk decays by a factor > 10 during 
the course of the simulation, the decay of Em proceeds much 
slow er, only by a fac tor of ~ 2. Similar to the incompressible 
case (lBiskampll2003h . in supersonic turbulence the energy ra- 
tio r = Ey^/Em is not constant but decreases with time. While 
Ms decreases by a factor of ~ 3 from 9 to 2.6, Ma shows 
a 2.5 X drop from ^ 4.5 to ^ 1.8. These differences in the 
decay rates of kinetic and magnetic energy as well as in the 
behavior of Ms and Ma can be understood as consequences 
of self-organization, i.e., the relaxation of the turbulence to- 
ward s an asyrnptotic static force-free minimum-energy state 
(e.g.. lBiskam"pll2003h . 

5.2. Small-scale Kinetic and Magnetic Diagnostics 

Another way to look at the effects of numerical dissipation 
is to analyze the time evolution of enstrophy r2=i<^|Vxup), 

" Note that, strictly speaking, the ILES approach involved here does not 
imply convergence as grid resolution improves since the effective Reynolds 
number is ultimately a function of the grid size (e.g., Kritsuk et al. 2006). At 
the same time, an asymptotic regime corresponding to Re = cxd can potentially 
be reached relatively early, at large, but still finite Reynolds numbers. This 
is what we probably observe as grid resolution approaches 2048^ in super- 
Alfvenic simulations. 
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Fig. 3. — Power spectra of the velocity (left panels) and magnetic energy (right panels) on a 512 grid for flow snapshots 1, 3, and 10 at f = 0.02, 0.06, and 
0.2 (top-to-bottom), respectively. The velocity spectra are compensated with k^l^ , while there is no compensation for the magnetic energy spectra. Note that the 
ordinate scale is not always the same for different time instances. 



dilatation A = (| V • up), and that of the mean squared current 
density 7^ = (|V x Bp). These so-called "small-scale" quan- 
tities show strong oscillations at the grid scale. Their spectra 
are dominated by the high wave numbers, and their PDFs have 
extended exponential tails (e.g. JPorter et al.ll2002l) . They also 
usually display a very slow (if any) convergence with grid res- 
olution in ILESs due to a strong dependence on Re and Rm. 
These measures are related to the total viscous a nd Ohmic dis- 
sipati on rates within the periodic domain (e.g., iKritsuk et alJ 
l2007h . For instance, in nonmagnetized compressible turbu- 



lence, which is in many respects similar to the super- Alfvenic 
case considered here, the mean dissipation rate of the specific 
kinetic energy can be expressed as ek = -(Re)"' (217 + 4/3 A) 
(e.g^ lPan et 10120091) .^° Since the global dissipation rates 
for iiK are very similar for all the methods considered here 
(see Figure[rii, the relative ranking of their effective Reynolds 
numbers is fully determined by the value of 20 + 4/3 A. We 
can thus use the dissipation rates plotted in the left panel of 

Strictly speaking, this is only vaHd for compressible Navier-Stokes tur- 
bulence assuming zero bulk viscosity, i.e., for ideal monoatomic gases. 
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Fig. 5. — Ratio of dilatational-to-solenoidal power in velocity spectra, x(^)> for 'he first (left panel) and the last (right panel) flow snapshots from the 512^ 
simulations at t[ = 0.02 and t[o = 0.2, respectively. 



Figure |2] to determine the relative standing of these meth- 
ods in terms of their Regff. The new STAGGER code shows 
an outstanding result during the first half of the evolution, 
f € [0,0.1], when its Reeff exceeds that of PPML by up to 
a factor of - 1 .5. RAMSES, PLUTO, LL-MHD, ENZO, KT- 
MHD, FLASH, and ZEUS follow STAGGER and PPML in 
the order of decreasing effective Reynolds number See Ta- 
ble 1 for numeric values of 257+4/3 A from different codes at 
t = 0.02. 

We employ the same approach to get an assessment of the 
relative standing of these numerical methods in terms of their 
effective magnetic Reynolds number, Rmgff. Figure |2] right 
panel shows the mean-squared current density, J^, as a func- 
tion of time. The current density is sensitive to both eM and 
Rmeff, since eM ^ -(Rm)~'7^. We expect qualitatively the 
same dependency here as for ek and Re^ff, but with, perhaps, 
different order of the methods. Note, however, that the dis- 
sipation rates and saturated levels of magnetic energy are not 
the same for different methods as can be seen in Figure[T] right 
panel. PPML shows the highest Rmeff, followed by FLASH, 
KT-MHD, RAMSES, PLUTO, ZEUS, ENZO, LL-MHD, and 
STAGGER. Note that the order in which methods follow each 
other in the right panel of Figure |2] is the same as in the 
plot in Figure [1] see also Table 1 . Thus, Rmgff and are 



positively correlated. We will explore this correlation further 
in Sections 15.31 and 15.41 where we analyze the power spec- 
tra of kinetic and magnetic energy and measure the effective 
bandwidth of the methods. 

5.3. Power Spectra 

Figure |3] shows power spectra of the velocity and magnetic 
energy at f = 0.02, 0.06, and 0.2. The spectra obtained at 
512^ demonstrate a very good agreement with each other up 
to logj()A:/A:,nin ~ 1.2 and slightly diverge at higher wavenum- 
bers. This means that numerical dissipation strongly affects 
the scales smaller or equal to ^ 16 grid cells in supersonic 
turbulence simulations with our best methods. 

The new STAGGER shows a very extended spectrum at f = 
0.02 with an asymptotic slope of -5/3 up to logiQ^/fc^in < 2. 
This slope is not preserved, however, for the whole duration of 
the simulation. By f = 0.2, when after many integration time 
steps the sonic Mach number drops to ~ 2.5, the spectrum 
looses some high-A: power and progressively bends down at 
logio/t/;tmin ^ 1.5 leaving behind PLUTO, RAMSES, PPML, 
and LL-MHD. A close inspection of the velocity spectra 
shows that numerical diffusion in ZEUS, FLASH, KT-MHD, 
and ENZO is somewhat stronger than in the rest of the grid- 
based codes and affects the velocities at lower wavenumbers. 
Besides the new STAGGER, RAMSES, PLUTO, PPML, and 
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LL-MHD are the least diffusive codes. The magnetic energy 
spectra display similar variations, but the ranking of methods 
is different. Here FLASH and PPML show very low magnetic 
diffusivity, while RAMSES, ZEUS, STAGGER, LL-MHD, 
and ENZO are more diffusive and KT-MHD stays in between. 

5.4. Effective Spectral Bandwidth 

In order to highlight variations in the power spectra ob- 
tained with different meth ods, we follow the proced ure devel- 
oped in iLele et al] (|2009^ and iJohnsen et alj ( 1201 Oh for com- 
pressible Navier-Stokes turbulence at moderate Mach num- 
bers. We have seen that PPML is one of the least diffusive 
methods here, so we declare the 1024^ PPML solution filtered 
with a low-pass Gaussian filter down to 256^* to be "exact" 
and call it the reference solution. We then plot power spectra 
compensated by the spectrum of this reference solution for 
the first snapshot at f = 0.02, see Figure |4] We set tolerance 
at a level of ±25% from the reference solution and define the 
spectral bandwidth of a method as the fraction of the Nyquist 
frequency where the compensated spectrum deviates by more 
than 25% from the reference solution. While this definition is 
rather arbitrary, it helps to establish a convenient quantitative 
measure to assert the performance of numerical methods for 
the turbulence decay test, see Table 1 for numeric values. 

The left panel of Figure |4] shows the compensated veloc- 
ity spectra. STAGGER, PPML, and RAMSES have the high- 
est spectral bandwidth in velocity. The second position is 
shared by LL-MHD, PLUTO, KT-MHD, and ENZO. ZEUS 
and FLASH show a similar velocity bandwidth. RAMSES 
and KT-MHD show a small bump at log^^^^k / k^^m G [0.6, 1.1] 
reminiscent of the bottleneck effect, while the STAGGER and 
PPML spectra decrease monotonically. 

The right panel of Figure |4] shows the compensated mag- 
netic energy spectra. The situation here is quite different. 
First, unlike the velocity spectra, the magnetic energy spec- 
tra start to bend down from the reference solution rather early. 
This is expected because of the slow convergence of with 
grid resolution, as we discussed earlier FLASH and PPML 
demonstrate the highest bandwidth, KT-MHD is in the mid- 
dle, PLUTO, ZEUS, LL-MHD, ENZO, and RAMSES are 
very similar to each other and show a somewhat higher mag- 
netic diffusivity. At 256^', the spectral bandwidth of our best 
MHD codes is ^ 0.3 for the velocities and ^ 0.2 for the mag- 
netic energy. If we were dealing with properly converged 
solutions ob tained from direct numerical simulations as in 
iJohnsen et a l. ( 2010), that would mean that numerical dissipa- 
tion strongly affects the wavenumbers down to at least (0.2- 
0.3)A:n. This is not however exactly the case here due to the 
adopted ILES approach, see footnote [T9] 

5.5. Dilatational versus Solenoidal Modes 

We have discussed above how the numerical methods dif- 
fer in their kinetic and magnetic diffusivity. This aspect 
plays an important role in simulations involving small-scale 
turbulent dynamo. There Pm serves as a control param- 
eter and the dynamo would on ly operate at Pm > Pmcrit 
dBrandenburg & Nordlundll201 lb . 

In this section, we look at how different methods treat di- 
latational and solenoidal components of the velocity field on 
small scales. We decompose the velocity fields into the po- 
tential (curl-free) and rotational (solenoidal) components, u = 
Ud + Us, using Helmholtz decomposition. We compute power 
spectra, P(ud,k) and P(Us,k), and define the dilatational-to- 
solenoidal ratio as x{k) = P{Ud,k)/P(Us,k). Peculiarities in 



the small-scale xik) ratio have a potential to affect various 
turbulence statistics (e.g., the density PDF) and limit (or 
even eliminate) the extent of the inertial range in simulations. 
These features cannot be captured by either the small-scale 
kinetic diagnostics discussed in Section 15.21 or by the power 
spectra discussed in Section 15.31 We present x(^) for snap- 
shots 1 and 10 (f = 0.02 and 0.2) in the left and right panels 
of Figure |5] respectively. Table 1 gives numeric values for 
the average dilatational-to-solenoidal ratio, x(^/^min > 100), 
at wavenumbers above lOOfcmin in the 5 12-' models at f = 0.2. 

First, note a very good agreement between all the methods 
in the early snapshot at low wavenumbers, log[Q^/A:min < 1 .3, 
with x(^/^min = 10) ~ 0.47. An overall ratio of 1 : 2 is ex- 
pected for super- Alfvenic turbulence at h igh Mach numbers 
(iKritsuk et a!]l2010l:lFederrath et al.ll2010h . As the turbulence 
decays, the sonic Mach number drops down to Ms ^ 2.7 by 
t = 0.2 and the ratio decreases to xi^/^min = 10) w 0.33, as 
expected. In the inertial range, xik) is k nown to be a slowly 
decre asing function of the wavenumber (iKritsuk et al.ll2007l 
1201 Ol) and this behavior is nicely captured by most of the 
codes. There are slight differences in the xik) levels between 
different methods at t = 0.2 with ZEUS and FLASH being 
slightly ahead of the other codes in damping the dilatational 
modes at high wavenumbers. Otherwise, the results from dif- 
ferent methods look very similar, although ENZO, KT-MHD, 
and STAGGER start to deviate somewhat from the rest of 
the codes at relatively low wavenumbers, log^Qk/k^i^ ~ 1.5. 
Also KT-MHD and ENZO produce unusually high x val- 
ues at the Nyquist wavenumber For instance, KT-MHD has 
xi^N) ~ 1-25 and 1.5 at f = 0.02 and 0.2, respectively. This 
indicates that perhaps some spurious compressible fluctua- 
tions are present at scales k > A:n/8 in simulations carried 
out with these codes. The small-scale oscillations of the 
KT-MHD code are likely to be caused by the way the CT 
scheme is implemented (R. Kissmann, 2010, private commu- 
nication). The observed compressible artifacts can probably 
be reduc ed to a large extent by using the CT approach pro- 
posed bv lLondrillo & del Zannal (|20()4|) . 

6. DISCUSSION 

One is tempted to try to sort these nine methods into some 
well-ordered set. This is, however, an impossible task, as no 
single solver consistently outperformed all others on all diag- 
nostics. In this discussion we will restrict our focus to dis- 
cussing kinematic and magnetic dissipation, as measured by 
the diagnostics presented here. This leaves out other poten- 
tially important di a gnost ics, such as the loop advection test of 
iGardiner & Stonel (120051) . We are also ignoring other salient 
features, such as computational cost (in memory and time to 
solution), ease or feasibility of extending the solver to differ- 
ent physical or numerical scenarios, or quality of documen- 
tation, all of which go into the selection of a code package. 
The final result is that all codes performed reasonably well 
on the task. There is no single silver bullet that determines 
the performance of a given solver; good quality results can be 
achieved through a variety of means, and dissipation can be 
introduced in a variety of places. 

All MHD algorithms used in this work are extensions of 
a previously established hydrodynamic algorithm. In gen- 
eral, five basic features determine the operation of a numer- 
ical scheme; base method (most prominently spatial order of 
accuracy), MHD extension, artificial viscosity, time integra- 
tion, and directional splitting. In this section, we will clas- 
sify each code based on these parameters and discuss trends 
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within each feature. In Section 16.11 we discuss the spatial 
order of accuracy, which seems to b e the dominant factor in 
determining performance. In Section we will discuss ar- 
tificial viscosity and source terms that behave as a viscosity. 
These features also seem to have considerable impact on the 
dissipation properties, as one would expect. In Section 16.31 
we will classify and discuss other properties of the base hy- 
drodynamical scheme. In Section|631 we will discuss the per- 
formance of various MHD extension methods. In Section |675l 
we will discuss three closely related schemes and discuss how 
seemingly small details can dramatically effect the results. In 
Section l676l we will discuss directional splitting, time integra- 
tion, and reconstructed variables. These seem to have a less 
dramatic impact on the overall performance, as measured by 
these diagnostics. See Table 2 for a summary of these solver 
configuration details. 

We refer th e read er to the excellen t books by iTorol (Il999l) . 
iLanevI (Il998l) . and iLeVequd (l2002h and the method papers 
cited in Section 3 for the details of each numerical scheme. 
We will not be expanding on any details except where neces- 
sary. 



6.1. Spatial Order of Accuracy 

High spatial order of accuracy seems to be the salient fea- 
ture of the least dissipative codes, though as there are many 
factors in each method that can improve or degrade perfor- 
mance. STAGGER has the highest spatial order, 6, and this 
is reflected most notably in Figure 2, left panel, where its ef- 
fective Reynolds number is significantly higher than the other 
methods, and Figure 3, left panel, where the inertial range 
of the power spectrum extends much further than the others. 
The third-order methods are PPML, PLUTO, KT-MHD, and 
the electric field construction of FLASH. These four meth- 
ods show the highest magnetic spectral bandwidth, and are 
the top performers in the effective magnetic Reynolds num- 
ber and magnetic power spectrum. However, other effects, 
most likely viscosity, keep these third-order methods from 
having the lowest dissipation among all statistics. The re- 
maining methods (ZEUS, RAMSES, LL-MHD, ENZO, and 
the base hydro scheme of FLASH) are second order spatially. 
These codes tend to show more dissipation over the third- 
order methods. 

There are two notable exceptions to this trend. The first can 
be seen in the spectral bandwidth plot, in which RAMSES 
(second order spatially) performs better than some third-order 
methods, though this may be due to other effects (see Section 
16.5b . The second exception can be seen in the top curve of 
the effective Reynolds number, corresponding to the STAG- 
GER method. The initial conditions were generated with an 
early version of STAGGER, but continued with a version that 
used conservative variables and different settings for the arti- 
ficial viscosity. As both methods are sixth order spatially, the 
increase in effective Reynolds number demonstrates that it is 
not spatial accuracy alone that determines dissipation prop- 
erties. This will be discussed further in Section 16.21 Note 
that here we refer only to the formal spatial order of accu- 
racy for the reconstruction or interpolation for each scheme. 
The actual convergence properties of each scheme, once time 
integration, spatial reconstruction, etc. have been taken into 
account, must be measured as a function of time and/or space 
resolution. This is beyond the scope of this work. 



6.2. Artificial Viscosity and Source Terms 

It is quite typical for numerical schemes to include some 
form of artificial viscos ity in order to avoid n ume rical instabil- 
ities. I n the case of the lPowell et al.l (Il999l) and lDedner et alJ 
(I2OO2I) MHD schemes, source terms proportional to V • B are 
included to constrain the effects of divergence, which while 
not the same kind of dissipation still have a dissipative ef- 
fect. In this suite of simulations, viscosity treatments can 
be broken coarsely into three categories: artificial viscosity, 
V • B motivated diffusivity, and exclusively numerical viscos- 
ity. ExpHcit terms are included in STAGGER, KT-MHD, 
ZEUS. Terms due to V • B treatments are included in PLUTO 
and ENZO. The four remaining codes have no explicit artifi- 
cial viscosity, and dissipation is only due to the scheme itself 
(these are FLASH, PPML, RAMSES, and LL-MHD). 

One naively expects that codes with explicit viscosity terms 
will have somewhat more dissipation than those without. 
However, this is only loosely seen in the results, and it is 
difficult to disentangle dissipative terms from other code dif- 
ferences. STAGGER gives the most noticeable example of 
the effects of dissipation, namely the large gap between its 
velocity dissipation, which is quite low, and its magnetic dis- 
sipation, which is quite a bit higher than other codes on most 
metrics. It is also possible that the fine tuning of the magnetic 
and kinematic artificial diffusivity, which has maximized the 
apparent inertial range, has altered the non-local coupling of 
MHD waves in a manner that still leaves the dissipation rel- 
atively high in the inertial range. It is reasonable to isolate 
codes based on spatial order of accuracy in order to com- 
pare viscosity results. Among the third-order codes, PPML, 
with no explicit viscosity, tends to have lower dissipation than 
PLUTO, which has V • B = motivated source terms, which in 
turn tends to have lower dissipation than KT-MHD. Then iso- 
lating the second order methods, the trend somewhat contin- 
ues, though less robustly. ENZO and ZEUS tend to show the 
most velocity dissipation, as measured by effective Reynolds 
number or velocity bandwidth. However, ENZO is the only 
second spatial order code with explicit magnetic dissipation, 
and it shows more power in the magnetic power spectrum than 
LL-MHD, which has none. RAMSES shows the lowest dis- 
sipation among the second order codes in all metrics except 
for magnetic spectral bandwidth, in which ZEUS is slightly 
higher 

One method (FLASH) includes terms proportional to the 
longitudinal derivative of the magnetic field. These terms are 
typically omitted from the derivation as they are identically 
zero in the one-dimensional version of the equations. 

6.3. Base Methods 

Eulerian hydro schemes fall broadly into two categories: 
finite-volume and finite difference. In loosest terms, finite 
volume schemes approximate the integral form of the con- 
servation law, while finite difference terms approximate the 
differential form of the conservation law. Three of the grid- 
based codes compared here are finite difference: STAGGER, 
KT-MHD, and ZEUS. The other six are finite-volume meth- 
ods (ENZO, FLASH, LL-MHD, PLUTO, PPML, and RAM- 
SES.) Between finite volume and finite-difference, there is no 
correlation with performance. This is best illustrated in the 
left panel of Figure 2. The code with the highest effective 
Reynolds number is STAGGER, and with the lowest is ZEUS, 
and both are finite difference methods. Here we will discuss 
some common features within each category of methods. 
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TABLE 2 

Solver Design Specieications eor the Eulerian Methods" 



Name 


Base Scheme'' 


Spatial Order*^ Source Terms'* 


MHD^ 


Time Integration* 


Directional Splitting^ 


ENZO 


FV, HLL 


2nd 


Dedner 


Dedner 


2nd-order RK 


Direct 


FLASH 


FY, HELD 


2nd 


1 1 Derivative 


3rd-order CT 


Forward Euler 


_L Reconstruction 


KT-MHD 


FD, CWENO 


3rd 


KT 


3rd-order CT 


4th-orderRK 


Direct 


LL-MHD 


FV, HELD 


2nd 


None 


Athena CT 


Forward Euler 


SpHt 


PLUTO 


FV, HELD 


3rd 


Powell 


Powell 


3rd-order RK 


Direct 


PPML 


FV, HEED 


3rd 


None 


Athena CT 


Forward Euler 


_L Reconstruction 


RAMSES 


FV, HEED 


2nd 


None 


2D HELD CT Forward Euler 


_L Reconstruction 


STAGGER FD, Stagger 


6th 


Tensor 


Staggered CT 3rd-order Hyman Direct 


ZEUS 


FD, van Leer 


2nd 


von Neumann 


MOC-CT 


Forward Euler 


Split 



" See Section|3]and the indicated sections on each topic for more information. 
Base method. FD for finite differ ence, FV for finite volume. FV techniques have the Riemann solver listed. Section |63] 
Spatial order of accuracy. Section |6T| 

Artificial Viscosity, Section "ii Derivative" indicates presence of terms proportional to the longitudinal derivative of the magnetic field. 
MHD method, Section lSH 



'Time integration method, Section l6.6.3l 
8 Multidimensional technique, Section |6.6.2| "± Reconstruction" indicates presence of transverse derivatives in the interface reconstruction. 



6.3.1. Finite -difference Methods 

One of the curses of numerical fluid dynamics is the bat- 
tle between accuracy and stability. These seem to be felt 
somewhat more strongly by the three finite-difference codes. 
ZEUS tends to be more dissipative than other methods, even 
though it is formally second order spatially (see, for instance, 
the effective Reynolds number in Figure 2, left, or the left col- 
umn of Figure 3). STAGGER has the highest effective fluid 
Reynolds number, but the lowest effective magnetic Reynolds 
number; we believe this to be a result of the tensor viscosity 
and its subtle relationship to the (scalar) magnetic diffusiv- 
ity. The KT-MHD method suffers from excessive small-scale 
compression, likely due to the fact that CWENO schemes 
are only essentially non-oscillatory, trading the possibility of 
small numerical oscillation near shocks for very high quality 
results in smooth regions. 

6.3.2. Finite-volume Methods 

The six finite-volume methods (ENZO, FLASH, LL-MHD, 
PLUTO, PPML, and RAMSES) are all some fom of higher- 
order extension of Godunov's method. These methods have 
the advantage that they capture shock structures, in principle, 
exactly. These methods can be broken into two parts; recon- 
struction to the interface, and the Riemann solver 

A wide array of Riemann solvers exist in the literature, but 
those used in this work are of two families. Roe and HLL. It 
is expected from other tests that Roe will perform the best, 
though it is subject to instabilities at low density and high 
Mach numbers, and HELD will perform the best of the HLL 
methods, as it captures more of the eigenstructure of the equa- 
tions than HLL. However, there does not seem to be a corre- 
lation between dissipation and choice of Riemann solver that 
cannot be explained by other mechanism. This does not say 
that there is no difference, merely not one that can be identi- 
fied by these data. 

The interface reconstruction techniques vary widely among 
the six schemes, and can primarily be characterized by details 
discussed in other sections. Namely, order of reconstruction, 
directional splitting, time integration, and explicit viscosity 
terms. They will not be discussed further here. 



6.4. MHD Methods 

Any MHD algorithm is essentially an established hydro- 
dynamic algorithm with modifications to include the Lorentz 
force in the momentum equation, the induction equation, and 
some treatment to minimize the divergence of the magnetic 
field. In all cases in this paper, the Lorentz force is incor- 
porated into the momentum equation directly (rather than 
through, say, vector potentials). Two of the codes (ENZO 
and PLUTO) use non-exact divergence preservation, namely 
both treat an extra wave, in V • B. These methods also include 
source terms to Equations (l)-(3) that are set to zero in most 
methods. This extra wave is advected and damped in ENZO, 
while it is simply advected with the fluid velocity in PLUTO 
and serves to mitigate singularities in the three-dimensional 
linearized Jacobian. The rest use a variant of the CT method, 
wherein the electric field and magnetic field are treated at the 
zone edge and face, respectively, which allows the solenoidal 
constraint to be kept zero to machine precision through the 
curl operator 

One naively expects the two approximate divergence meth- 
ods to have somewhat higher dissipation than other codes, as 
the primary driver is dissipation. This is to some extent seen in 
the data, though PLUTO does not suffer much from this effect 
as it still has quite high fluid and magnetic Reynolds numbers. 
ENZO, on the other hand, seems to have more dissipation, and 
based on its similarity to other spatially second-order codes, 
the V • B wave seems to be a likely culprit. 

Among the CT-based schemes, the results seem to be 
dominated by first spatial order then on reconstructed vari- 
able. PPML and LL-MHD both use the electric field re- 
construction technique d escribed in the Athena method of 
iGardiner & Stonel (l2005l) . but PPML is spatially third-order, 
so has a higher magnetic Reynolds number It should be 
noted that only the electric field reconstruction of the Athena 
method is used by these methods. FLASH also uses third- 
order reconstruction, and also has an extremely large mag- 
netic spectral bandwidth. LL-MHD, ZEUS, and RAMSES, 
on the other hand, are all spatially second-order, but ZEUS 
uses MOC, which uses the characteristic fields, and RAM- 
SES solves a second Riemann problem; both of which prove 
to better capture the electric field than the primitive variable 
reconstruction used in LL-MHD. 
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6.5. Three Closely Related Codes 

An interesting subset of codes to examine are LL-MHD, 
RAMSES, and FLASH. These three codes are the most sim- 
ilar in terms of their components, and serve to illustrate how 
small differences in method details can cause significant dif- 
ferences in the performance. Each of the three codes uses 
the second-order MUSCL-Hancock reconstruction-evolution 
scheme for computation of interface states, the HELD Rie- 
mann solver, forward Euler time integration, and a higher- 
order CT method. Given all the similarities, the differences 
in performance of the three codes are surprising. This is best 
shown in the two spectral bandwidth plots in Figure 4. The 
magnetic bandwidth of FLASH, in Figure 4, right, is the high- 
est of all available codes, as measured by the wavenumber 
at which the spectrum crosses 75%, with log[QA:/A;niin=1.4. 
RAMSES and LL-MHD are significantly lower, both with 
log[QA:/A;niin=l.l. The spatial order of accuracy in the elec- 
tric field computation is a clear culprit. FLASH uses a third- 
order central-difference reconstruction of the electric fields 
from the Riemann solver RAMSES and LL-MHD, on the 
other hand, both use spatially second-order methods, with 
RAMSES using a novel two-dimensional Riemann solver, 
and LL-MHD using the Athena method. This shows the 
importance of spatial reconstruction in capturing flow fea- 
tures. The velocity bandwidth, in Figure 4, left, is a com- 
pletely different story: RAMSES is at the high end of the 
codes, with logjQA:/A:min=1.5, LL-MHD is in the middle, with 
log[oA:/A:min=1.4, but FLASH is the lowest of the Eulerian 
codes, with log[o^/^min=l-l- This is most interesting, since 
the base solvers for each of the codes are nearly identical. The 
biggest difference here is the treatment of directional splitting. 
RAMSES and FLASH are both directionally unsplit, incor- 
porating transverse derivatives of the Jacobian in the interface 
reconstruction as discussed in Section 16.6.21 while LL-MHD 
is split using Strang splitting, and does not include transverse 
derivatives. This alone does not explain the ordering, as LL- 
MHD and RAMSES perform quite similarly in many veloc- 
ity statistics. The only other major algorithmic difference is 
the inclusion of the longitudinal magnetic derivatives in the 
FLASH interface reconstruction. It is not obvious that these 
terms would cause diffusion in the manner observed, though 
they will affect the reconstruction of the interface states. Fi- 
nally, each of these methods (indeed all methods described 
here) include a number of nonlinear switches that determine 
behavior near shocks, among other things, that have not been 
explicitly described. Further investigation is required to iso- 
late these finer details. 

Due to the tight coupling between the velocity and the 
magnetic field in both the momentum and induction equa- 
tions, it would not be surprising for the velocity and magnetic 
statistics to be coupled, even perhaps in an inverse manner, 
through either energy conservation or mode coupling. Thus 
the higher spatial order used in the FLASH magnetic recon- 
struction may, for instance, be more efficient at transferring 
kinetic energy to magnetic energy. More study is required to 
definitively pinpoint the cause of differences between these 
three codes, but it illustrates the effect of seemingly minor 
details having substantial results in the behavior of a code. 

An additional point of interest in the RAMSES behavior 
is the excess power seen in the spectral bandwidth plot at 
logjo^/fcrnin G [0.8, 1.2]. This seems to be a manifestation of 
what in pure hydrodynamic turbulence is referred to as the 
bottleneck. This is typically not seen in simulations of MHD 



turbulence at a 512 resolution, presumably due to additional 
effects of non-local MHD mode coupling that allows energy 
to be more efficiently transferred to smaller scales. As RAM- 
SES has a relatively low Prandtl number, it is possible that 
this extra energy transfer is not as efficient as in other codes, 
causing somewhat inflated spectral bandwidth. 

6.6. Other Solver Details 

There are several other solver design specifications that 
have received considerable attention over the years. Here we 
present a discussion of some of the major solver options that 
have been examined. While each may be crucial in its own 
right, are not dominant factors determining the dissipation 
properties studied here. 

6.6.1. Evolved Variables 

MHD can be described by three distinct sets of vari- 
ables; primitive variables (/9,u,B,/7tot); conserved variables 
(PjpUjBjEtot), and the characteristic variables, which are 
the eigenvectors of the Jacobian of the equations, and in some 
ways the most physically relevant form of the variables. It has 
been shown that in some cases working with the cha racteris- 
tic va riables gives superior results to the other two (iBalsaral 
l2004l) . There is some evidence that bears this out in these data. 
For instance, the magnetic behavior of ZEUS, in which MOC 
traces characteristics to compute the electric field, is generally 
less dissipative than LL-MHD, which uses primarily primitive 
variables. However, this is not universally the case, and other 
factors may prove more important. Such is the case in ve- 
locity performance of FLASH vsersus RAMSES, which use 
spatial limiting on characteristic and primitive variables, re- 
spectively. 

6.6.2. Directional Splitting 

Computational algorithms have a long history of being de- 
veloped as one-dimensional methods. They then must be ex- 
tended by some manner to three dimensions. There are es- 
sentially three categories of multidimensional techniques em- 
ployed by the codes in this study; directly unsplit, direction- 
ally split, and transverse flux methods. 

The two directionally split methods (LL-MHD and ZEUS) 
employ sequential one-dimensional solves along each coordi- 
nate axis, wherein the partial update of one sweep is used as 
the initial data for the following sweep. The order of sweeps 
in both methods is permuted to reduce error In both meth- 
ods, the electric field is computed after the three sweeps are 
finished. 

The four "directly unsplit" methods are STAGGER, KT- 
MHD, PLUTO, and ENZO. The first two do not rely on 
strictly one-dimensional techniques, so they employ fully 
three-dimensional evolution by repeated application of the in- 
terpolation and derivative operators. The ENZO and PLUTO 
methods use the Godunov method, which is strictly speak- 
ing one-dimensional as will be discussed below. It applies 
the algorithm in an unsplit fashion, with the initial state for 
Riemann solves coming from the same data for all three di- 
mensions. It incorporates multidimensional properties of the 
flow by way of the Runge-Kutta integration. 

Another unsplit technique, dubbed "transverse reconstruc- 
tion" here, is used to incorporate three-dimensional terms into 
the finite-volume methods. Godunov's method is, strictly 
speaking, one-dimensional, and does not lend itself di- 
rectly to multidimensional techniques. The underlying one- 
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dimensional method follows three basic steps; reconstruc- 
tion of two interface states at each zone boundary, followed 
by solution of the Riemann problem at the zone bound- 
ary, and finally using that solution to compute and differ- 
ence fluxes at the interface to update the field. Three of 
the schemes (FLASH, PPML, and RAMSES) introduce the 
multidimensional terms in the reconstruction of the interface 
state, through the addition of terms approximating gradients 
of the transverse fluxes, dFy/dy. The techniques vary slightly 
between the four methods. FLASH and RAMSES use a lin- 
earization of the transverse flux gradient, AydU/dy to com- 
pute the half step advance in time. Here, Ay is the Jacobian 
of the flux, and the derivative is approximated with mono- 
tonized central differences. PLUTO uses a full reconstruction 
and Riemann solve in the transverse direction. PPML also in- 
cludes a linearization of the transverse flux, though it is incor- 
porated slightly differently, with the transverse flux gradient 
introduced in the characteristic invariants, and uses character- 
istic direction filtering for upwinding the derivative. 

While it is often believed that directional sweeping is a 
detriment to the solution quality, the metrics presented in this 
work do not show a clear correlation between the different 
multidimensional techniques. 

6.6.3. Time Integration 

In principle, the order of the spatial and temporal integra- 
tion should be the same, otherwise the convergence properties 
of the scheme will be reduced to the lower of the two. How- 
ever, time integration in these cases seems to be swamped by 
other effects. ENZO is of higher order in time than RAMSES, 
but significantly more dissipative. Similarly, PLUTO is higher 
order in time than PPML, but typically has higher dissipation, 
as well. 

7. SUMMARY AND CONCLUSIONS 

We have compared nine numerical MHD codes on a decay- 
ing supersonic, super- Alfvenic turbulence test problem with 
conditions similar to star-forming molecular clouds in the 
Galaxy. The codes ENZO, FLASH, KT-MHD, LL-MHD, 
PLUTO, PPML, RAMSES, STAGGER, and ZEUS, described 
in detail in Section 2, employ a variety of numerical algo- 
rithms of varying order of accuracy, multidimensional and 
time integration schemes, shock capturing techniques, and 
treatment of the solenoidal constraint on the magnetic field. 
Together, they represent a majority of the MHD codes in use 
in numerical astrophysics today and therefore sample the cur- 
rent state of the art. The work presented in this paper is the 
largest, most comprehensive MHD code comparison on an 
application-like test problem to date. 

The codes were compared using both integrated and spec- 
tral measures of the velocity and magnetic fields. All nine 
Eulerian codes agreed surprisingly well on the kinetic en- 
ergy decay rate (Figure 1, top left), which indicates both 
the rob ustness of published predictions ("Mac Low et al.lfT998l : 
IStone e t al. 1998; Lemaster & Stone 2009) as well as the in- 
adequacy of this particular metric as a discriminant among 
methods. All nine Eulerian codes likewise agreed on the mag- 
netic energy decay rate (Figure 1, top right), but varied on the 
amplitude of the peak magnetic energy as this proved sensitive 
to the effective magnetic Reynolds number of the simulation, 
which depends on the numerical dissipation of the method. 

To move beyond simple global energy diagnostics, small 
scale kinetic and magnetic field diagnostics were introduced 



in order to empirically measure the effective fluid and mag- 
netic Reynolds numbers of the various codes. These diag- 
nostics are based on analytically motivated combinations of 
the volume integrated fluid enstrophy, dilatation, and square 
of the electric current density (Figure 2). They proved more 
revealing about the numerical dissipation present in the var- 
ious methods, and motivated a closer investigation using 
power spectra of the velocity and magnetic fields. Regarding 
the latter, the concept of effective spectral bandwidth (ESB) 
was introduced as a quantitative metric for code comparison. 
The effective spectral bandwidth is defined as the width in 
wavenumber space where the numerical results do not deviate 
from a reference solution (typically, a higher resolution sim- 
ulation) by more than 25%. The ESB was measured for both 
the velocity and magnetic power spectra for all nine codes at 
reference times during the decay. A detailed comparison of 
ESB's leads to several general conclusions and observations. 

1 . All codes gave qualitatively the same results, implying 
that they are all performing reasonably well and are use- 
ful for scientific investigations. 

2. No single code outperformed all the others against 
all metrics, although higher-order-accurate methods do 
better than lower-order-accurate methods in general. 
The lack of a clear winner stems from the fact that a 
single MHD code is a combination of many different al- 
gorithms representing specific design choices, and that 
many combinations are possible. 

3. The spatial order of accuracy is the primary determinant 
of velocity spectral bandwidth and effective Reynolds 
number Higher spatial order correlates to higher spec- 
tral bandwidth. The sixth-order code STAGGER is su- 
perior to the third-order codes PPML, PLUTO, KT- 
MHD and FLASH, which are superior to the second- 
order codes ZEUS, LL-MHD, and ENZO. 

4. Codes with high velocity spectral bandwidth do not 
necessarily have high magnetic spectral bandwidth. For 
example, the STACiGER code has the highest velocity 
ESB but the lowest magnetic ESB. The magnetic ESB 
is sensitive to the spatial order of accuracy of the elec- 
tric field computation, and is higher in methods that in- 
terpolate on characteristic variables as opposed to prim- 
itive variables. 

5. The use of explicit artificial viscosity to stabilize shock 
waves reduces the velocity spectral bandwidth relative 
to methods that do not use artificial viscosity, such as 
Godunov methods. 

6. The use of explicit divergence cleaning reduces the 
magnetic spectral bandwidth relative to codes that pre- 
serve the solenoidal condition on B exactly (CT meth- 
ods). 

7. Other algorithmic choices such as finite-difference ver- 
sus finite-volume discretization, directionally split ver- 
sus unsplit updates of the conservations laws, and order 
of accuracy of the time integration are less well cor- 
related with the performance metrics, and therefore ap- 
pear to be less important in predicting a code's behavior 
on MHD turbulence. 

Observations about specific codes are as follows: 
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• The best performers overall are PPML, FLASH, 
PLUTO, and RAMSES based on velocity and magnetic 
Reynolds numbers and spectral bandwidths. 

• The highest fluid Reynolds number was obtained with 
the STAGGER code. 

• The highest magnetic Prandtl number was obtained 
with the FLASH code. 

• FLASH is somewhat more diffusive on the hydro part 
than its magnetic part, and the reverse is true for the 
RAMSES code. 

• The dilatation velocity power spectra of KT-MHD and 
ENZO exhibit problematic behavior on small scales 
that is likely related to the ways these codes maintain 
V-B = 0. 

The best performing codes employ a consistently high or- 
der of accuracy for spatial reconstruction of the evolved fields, 
transverse gradient interpolation, conservation law update 
step, and Lorentz force computation. Three of the four em- 
ploy divergence-free evolution of the magnetic field using the 
CT method, and all use little to no explicit artificial viscosity. 
These would seem to be guidelines for the development of fu- 
ture schemes. Codes which fall short in one or more of these 
areas are still useful, but they must compensate higher numer- 
ical dissipation with higher numerical resolution. A new class 
of nearly Lagrangian methods for hydrodynamics has recently 
emerged which uses a moving mesh based on Voronoi cells 
([Springel 2010). It remains to be seen if this approach can be 
generalized to MHD while retaining the beneficial elements 
of successful Eulerian schemes. 
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